---
title: "4 - H4"
author: "Arturo Bertero"
date: ""
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
#clean
rm(list = ls())

# Libraries
library(pacman)
p_load(tidyverse, here, EGAnet, psychTools, conflicted, haven, lavaan, qgraph, forcats,
       stargazer, mgm, corclass, ggcorrplot, ggpubr, nnet, ggeffects, magick, semTools,
       sjPlot, patchwork) 

#notation
options(scipen=999)

#Conflicts
conflicts_prefer(dplyr::select)
conflicts_prefer(dplyr::filter)
```

# Input

```{r}
#Load data
IPBS = readRDS((here("Input", "IPBS_with_cca.rds")))

#Filter smaller dataset
att = IPBS %>% 
  select(c(abort:redis, globa, immig)) 

```

# Processing

```{r}
# Define shortnames, longnames, and groupings 
shortnames <- names(att)
longnames <- c("Abortion","Euthanasia","Same-sex marriage",
               "Redistribution", "Globalization","Immigration")

# Obj for dimensions
Totalgroup_comm <- list(
  "Cultural" = c(1:3),
  "Economic" = c(4:6)
)

Totalgroup_cols <- c("#966FD6B3","#008B8BB3")  

#remove labels
# att = sapply(att, haven::zap_labels)
```


## Multinomial regression

```{r}
# Drop NA and unused levels
data_m <- IPBS %>%
  filter(!is.na(cca), !is.na(vote_cat))

# Drop unused cca levels and relabel them as 1, 2, 3
data_m$cca <- factor(data_m$cca)
data_m$cca <- fct_drop(data_m$cca)
levels(data_m$cca) <- c("1", "2", "3")

# Ensure i have factors
data_m$vote_cat <- factor(data_m$vote_cat)
data_m$educ_cat <- factor(data_m$educ_cat)
data_m$pol_int_cat <- factor(data_m$pol_int_cat)
data_m$hh_income_cat <- factor(data_m$hh_income_cat)
data_m$sex <- factor(data_m$sex)
data_m$L_R_cat <- factor(data_m$L_R_cat)
```

### M1: pol int on cca

```{r}
# M1: pol int on membership
m1 = multinom(cca ~ pol_int + educ_cat + hh_income + sex + age, data = data_m)

#Tab
tab_model(m1, p.style = "stars",
          title = "", show.intercept = FALSE,
          show.loglik = TRUE, show.aic = TRUE, 
          file = here("Output", "Supplement", "Tab_6.doc"))
```

```{r}
# Get predicted probabilities 
pred_m1 <- ggeffect(m1, terms = c("pol_int"))

#Modify values 
pred_m1$response.level = c("Class 1\nIdeologues", "Class 1\nIdeologues", "Class 1\nIdeologues", "Class 1\nIdeologues",
                           "Class 2\nAlternatives", "Class 2\nAlternatives", "Class 2\nAlternatives", "Class 2\nAlternatives", 
                           "Class 3\nFragmented", "Class 3\nFragmented", "Class 3\nFragmented", "Class 3\nFragmented")

#Plot
plot_pint_by_cca <- ggplot(pred_m1, aes(x = x, y = predicted, group = response.level,
                                     color = response.level, fill = response.level)) +
  geom_line(size = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.15, color = NA) +
  scale_color_brewer(palette = "Dark2", name = "Belief Class") +
  scale_fill_brewer(palette = "Dark2", name = "Belief Class") +
  labs(
    x = "Political Interest",
    y = "Predicted Probability",
    title = ""
  ) +
  facet_wrap(~response.level, nrow = 1) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.title = element_text(face = "bold"),
    axis.text.x = element_text(hjust = 1),
    legend.position = ""
  )

#Save
ggsave(here("Output", "Article", "Figure_5.jpeg"),
       plot_pint_by_cca,
       dpi = 600,
       height = 8, width = 12)

#Save as tiff for publication
ggsave(
  here("Output", "Article", "Figure_5_tiff.tiff"),
  plot_pint_by_cca,
  dpi = 600,                 
  width = 12,
  height = 8,
  device = ragg::agg_tiff,   
  compression = "lzw",       
  bg = "white"               
)
```

### M2: L-R cat on cca

```{r}
# M2: sociodem + L-R on membership
m2 = multinom(cca ~ L_R_cat + pol_int + educ_cat + hh_income + sex + age, data = data_m)

#Tab
tab_model(m2, p.style = "stars",
          title = "", show.intercept = FALSE,
          show.loglik = TRUE, show.aic = TRUE, 
          file = here("Output", "Supplement", "Tab_7.doc"))
```

```{r}
# Get predicted probabilities for lr by CCA
pred_m2 <- ggeffect(m2, terms = "L_R_cat")

pred_m2$response.level = c("Class 1\nIdeologues", "Class 1\nIdeologues", "Class 1\nIdeologues", "Class 1\nIdeologues",
                           "Class 2\nAlternatives", "Class 2\nAlternatives", "Class 2\nAlternatives", "Class 2\nAlternatives", 
                           "Class 3\nFragmented", "Class 3\nFragmented", "Class 3\nFragmented", "Class 3\nFragmented")

# Ensure L_R_cat is ordered logically
pred_m2$x <- factor(pred_m2$x, levels = c("Left", "Center", "Right", "No lr"))


# Plot
plot_cca_by_lr <- ggplot(pred_m2, aes(x = x, y = predicted, group = response.level,
                                     color = response.level, fill = response.level)) +
  geom_point(position = position_dodge(width = 0.5), size = 2) +
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high),
    position = position_dodge(width = 0.5),
    width = 0.2
  ) +
  scale_color_brewer(palette = "Dark2", name = "Belief Class") +
  scale_fill_brewer(palette = "Dark2", name = "Belief Class") +
  labs(
    x = "Left-Right",
    y = "Predicted Probability",
    title = "M2: Predicted Probability of Belief System Membership by L-R"
  ) +
  facet_wrap(~x, nrow = 1) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.title = element_text(face = "bold"),
    axis.text.x = element_text(hjust = 1),
    legend.position = ""
  )

```


### M3 vote on cca

```{r}
# M3: sociodem + vote on membership
m3 = multinom(cca ~ vote_cat + pol_int + educ_cat + hh_income + sex + age, data = data_m)

#Tab
tab_model(m3, p.style = "stars",
          title = "", show.intercept = FALSE,
          show.loglik = TRUE, show.aic = TRUE, 
          file = here("Output", "Supplement", "Tab_8.doc"))
```

```{r}
# Get predicted probabilities for lr by CCA
pred_m3 <- ggeffect(m3, terms = "vote_cat")

pred_m3$response.level = c("Class 1\nIdeologues", "Class 1\nIdeologues", 
                           "Class 1\nIdeologues", "Class 1\nIdeologues",
                           "Class 1\nIdeologues", "Class 2\nAlternatives",
                           "Class 2\nAlternatives", "Class 2\nAlternatives", 
                           "Class 2\nAlternatives","Class 2\nAlternatives", 
                           "Class 3\nFragmented", "Class 3\nFragmented", 
                           "Class 3\nFragmented", "Class 3\nFragmented",
                           "Class 3\nFragmented")

pred_m3$x = c("Right", "M5S", "No Vote","Other", "Left",
                           "Right", "M5S", "No Vote","Other", "Left",
                          "Right", "M5S", "No Vote","Other", "Left")


# Ensure L_R_cat is ordered logically
pred_m3$x <- factor(pred_m3$x, levels = c("Left", "M5S", "Other", "Right", "No Vote"))

# Plot
plot_M3 <- ggplot(pred_m3, aes(x = x, y = predicted, group = response.level,
                                     color = response.level, fill = response.level)) +
  geom_point(position = position_dodge(width = 0.5), size = 2) +
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high),
    position = position_dodge(width = 0.5),
    width = 0.2
  ) +
  scale_color_brewer(palette = "Dark2", name = "Belief Class") +
  scale_fill_brewer(palette = "Dark2", name = "Belief Class") +
  labs(
    x = "Vote Choice",
    y = "Predicted Probability",
    title = "M3: Predicted Probability of Belief System Membership by Vote Choice"
  ) +
  facet_wrap(~x, nrow = 1) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.title = element_text(face = "bold"),
    axis.text.x = element_text(hjust = 1),
    legend.position = "bottom"
  )


#Save
ggsave(here("Output", "Supplement", "Fig_3.jpeg"),
       plot_M3,
       dpi = 600,
       height = 10, width = 12)
```


### Assemble

```{r}
#Together
M2_M3 = plot_cca_by_lr / plot_M3

#Save
ggsave(here("Output", "Supplement", "Fig_3.jpeg"),
       M2_M3,
       dpi = 600,
       height = 14, width = 12)
```


## M4: is cca membership more influenced by vote or by lr?

```{r}
# Model 4: both L_R_cat and vote_cat
m4 <- multinom(cca ~ L_R_cat + vote_cat + pol_int + educ_cat + hh_income + sex + age, data = data_m)

# Compare AICs
aic_comp = AIC(m1, m2, m3, m4)
bic_comp = BIC(m1, m2, m3, m4)

# reg tab
tab_model(m4, p.style = "stars",
          title = "",
          show.intercept = FALSE,
          show.loglik = TRUE,
          show.aic = TRUE,
          file = here("Output", "Supplement", "Tab_9.doc"))

```







